{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scripts.forest import Forest\n",
    "from sklearn.ensemble import IsolationForest\n",
    "from sklearn import mixture\n",
    "from scipy.linalg import qr\n",
    "import scipy.stats as scp\n",
    "from math import cos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#n_pts per Gaussian\n",
    "def sample_points(n_dimensions, n_pts, n_zeros):\n",
    "    n_ones = n_pts - n_zeros\n",
    "    pts = np.zeros((n_dimensions,n_pts))\n",
    "    for i in range(n_dimensions):\n",
    "        pts[i,:n_ones] = (2*np.random.randint(2, size=n_ones))-1\n",
    "    return pts\n",
    "        \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_dimensions = 10\n",
    "n_pts = 1000\n",
    "n_zeros = 30\n",
    "pts = sample_points(n_dimensions, n_pts, n_zeros)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_iso(pts, sample_num):\n",
    "    rng = np.random.RandomState()\n",
    "    dim,n_p = np.shape(pts)\n",
    "    clf = IsolationForest(max_samples = sample_num, random_state = rng, contamination = 0.05, n_estimators= 20, behaviour = \"new\")\n",
    "    clf.fit(np.transpose(pts))\n",
    "    Y = clf.predict(np.transpose(pts))\n",
    "    iso_indices = []\n",
    "    for i in range(len(Y)):\n",
    "        if Y[i] == -1:\n",
    "            iso_indices.append(i)\n",
    "    return iso_indices\n",
    "\n",
    "def get_density(pts, sample_num):\n",
    "    kwargs = {'max_depth': 20, 'n_trees':20,  'max_samples': sample_num, 'max_buckets': 2, 'epsilon': 0.1, 'sample_axis': 1, 'threshold': 0}\n",
    "    forest = Forest(**kwargs)\n",
    "    forest.fit(pts)\n",
    "    gsw_indices, outliers, scores , pst = forest.predict(pts, 0.05)\n",
    "    return gsw_indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[30. 30. 30. 30. 30. 30.  0. 30.  0.  0. 30.  0.  0.  0.  0.  0.  0.  0.\n",
      "  0.]\n",
      "[60. 60. 60. 60. 60. 30. 30. 30.  0.  0. 30.  0.  0.  0.  0.  0.  0.  0.\n",
      "  0.]\n",
      "[90. 90. 90. 90. 90. 30. 30. 30.  0. 30. 30.  0.  0. 30.  0.  0.  0.  0.\n",
      "  0.]\n",
      "[120. 120. 120. 120.  90.  60.  30.  60.  30.  30.  60.   0.   0.  30.\n",
      "   0.   0.   0.   0.   0.]\n",
      "[150. 150. 150. 150. 120.  60.  30.  60.  30.  30.  60.   0.   0.  30.\n",
      "   0.   0.   0.   0.   0.]\n",
      "[180. 180. 180. 180. 150.  90.  60.  60.  60.  30.  60.   0.   0.  30.\n",
      "   0.   0.   0.   0.   0.]\n",
      "[210. 210. 210. 210. 150. 120.  90.  60.  60.  30.  60.   0.   0.  30.\n",
      "   0.   0.   0.   0.   0.]\n",
      "[240. 240. 240. 240. 180. 150.  90.  90.  60.  30.  60.   0.   0.  30.\n",
      "   0.   0.   0.   0.   0.]\n",
      "[270. 270. 270. 270. 210. 150.  90.  90.  60.  30.  90.   0.   0.  30.\n",
      "   0.   0.   0.   0.   0.]\n",
      "[300. 300. 300. 300. 240. 180. 120. 120.  60.  30.  90.  30.   0.  30.\n",
      "   0.   0.   0.   0.   0.]\n",
      "[30. 30. 30. 30. 24. 18. 12. 12.  6.  3.  9.  3.  0.  3.  0.  0.  0.  0.\n",
      "  0.]\n"
     ]
    }
   ],
   "source": [
    "n_ex = 19\n",
    "n_rep = 10\n",
    "res_iso = np.zeros(19)\n",
    "truth = range(970,1000)\n",
    "for j in range(n_rep):\n",
    "    for i in range(19):\n",
    "        iso_indexes = get_iso(pts,100 + 50*i)\n",
    "        res_iso[i] += len(set(iso_indexes).intersection(set(truth)))\n",
    "    print(res_iso)\n",
    "res_iso = res_iso / n_rep\n",
    "print(res_iso)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "30.0\n",
      "30.0\n",
      "30.0\n",
      "30.0\n",
      "30.0\n",
      "30.0\n",
      "30.0\n",
      "30.0\n",
      "30.0\n",
      "30.0\n",
      "30.0\n",
      "30.0\n",
      "30.0\n",
      "30.0\n",
      "30.0\n",
      "30.0\n",
      "30.0\n",
      "30.0\n",
      "30.0\n",
      "[30. 30. 30. 30. 30. 30. 30. 30. 30. 30. 30. 30. 30. 30. 30. 30. 30. 30.\n",
      " 30.]\n"
     ]
    }
   ],
   "source": [
    "res = np.zeros(19)\n",
    "truth = range(970,1000)\n",
    "for i in range(19):\n",
    "    gsw_indexes = get_density(pts,100 + 50*i)\n",
    "    res[i] = len(set(gsw_indexes).intersection(set(truth)))\n",
    "    print(res[i])\n",
    "print(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xd8VFX6+PHPk0JCSAFCEoEAAQkIKjU0EVGxICpNEBAF7K66oLsW1nVddr+ua9mf7lrWsgqoKEURRMWGCqIUCVWqgAQIUhIgIYQkpJzfH/fOMISETJJpSZ736zWvzNxy7jMl88y559xzxBiDUkopBRDk7wCUUkoFDk0KSimlnDQpKKWUctKkoJRSykmTglJKKSdNCkoppZw0KaiAIiL1ReQTEckWkQ/sZU+KSKaIHPB3fGURkb4isl1EjovIUD/FcKmIpJezrp+IbPN1TKpm0qSg/EJEFovIUREJK7VqBJAAxBpjRopIS+CPQEdjzDnVOF65X5oe8HfgZWNMpDFmvpeOUWXGmKXGmPb+jkPVDJoUlM+JSBLQDzDA4FKrWwG/GGOK7MctgcPGmEM+C7DyWgGb/B2EUp6gSUH5wzhgBTAdGO9YKCJ/A54ARtmnYu4Gvgaa2Y+n29v1FpFlIpIlIutF5FKXMhqLyDQR+c2uicwXkQbA5y7lHBeRZiLSU0RSReSYiBwUkefLC1hE7hSRHSJyREQWiEgze/lOoA3wiV1u6ZoPIpImIg+LyAYRyRWRt0QkQUQ+F5EcEVkkIo1ctv9ARA7Yp9C+F5HzXdYNEpHN9n77ROShcuKdaG+XWLqWZMfzkB1PtojMFpFwl/WPiMh++zW8Q0SMiLQt77VRtYwxRm968+kN2AHcC3QHCoEEl3VTgBkujy8F0l0eNwcOA4OwftRcaT+Os9d/BswGGgGhQP+yyrGXLQduse9HAr3LifdyIBPoBoQBLwHfu6xPA644y/NNw0qCCXb8h4A1QFcgHPgW+KvL9rcBUfax/g2sc1m3H+hn328EdCv9/LAS6xqX16T0a5gG/AQ0AxoDW4B77HUDgQPA+UAEMAOrRtfW358bvfnmFnLWjKGUh4nIxVinW+YYYzLtX9o3AS+4WcTNwEJjzEL78dcikgoMEpGvgGuw2iOO2uuXnKWsQqCtiDQxxmRifXGXZSww1Rizxn4OfwKOikiSMSbNzbhfMsYctPdfChwyxqy1H88DBjg2NMZMddwXkSn2sWKMMdl2zB1FZL39HI+6HEPs2k5P4DJ7+/K8aIz5zd7pE6CLvfxGYJoxZpPL8ce6+RxVLaCnj5SvjQe+sr+EAd7H5RSSG1oBI+1TR1kikgVcDDQFWgBHXBJCRW4H2gFbRWSViFxXznbNgN2OB8aY41i1k+aViPugy/28Mh5HAohIsIg8LSI7ReQY1q96gCb23xuwakm7RWSJiPRxKachcBfwzwoSAli1AYcTjuNjPde9Lutc76s6QGsKymdEpD7WL9Fgl+6lYUBDEelsjFnvRjF7gXeNMXeWUX5ToLGINDTGZJVafcZwwMaY7cAYEQkChgMfikisMSa31Ka/YSUjx3EaALHAPjfiraybgCHAFVgJIQarNiB2zKuAISISCtwPzMFKhtjb3QzMEZFhxpgfq3D8/UCiy+MW5W2oaietKShfGgoUAx2xTld0AToAS7Ean90xA7heRK62f1WH2w2picaY/VgNyv8VkUYiEioil9j7HQRiRSTGUZCI3CwiccaYEsCRRErKOOZM4FYR6WI3JD8FrKzEqaPKiAIKsGoiEfaxHPHWE5Gx9qmkQuBY6XiNMYuxTvd8JCI9q3D8OVjPtYOIRAB/qdrTUDWVJgXlS+OxzlfvMcYccNyAl4GxIlJhzdUYsxfrl/RjQAZWzeFhTn2Wb8E6774Vq0H3AXu/rVhf7r/ap52aYTWqbhKR48B/gNHGmLwyjrkI68txLtYv6XOB0VV8DSryDtapqn3AZs5s57gFSLNPLd1DGef7jTFfYzVWfyIi3SpzcGPM58CLwHdYHQIcxy+oTDmq5hJjdJIdpVTZRKQDsBEIM6euHVG1mNYUlFKnEZFhIhJmXzvxDPCJJoS6Q5OCUqq0u7FOve3EagP6nX/DUb6kp4+UUko5aU1BKaWUU427TqFJkyYmKSnJ32EopVSNsnr16kxjTFxF29W4pJCUlERqaqq/w1BKqRpFRHZXvJWePlJKKeVCk4JSSiknTQpKKaWcNCkopZRy0qSglFLKyWu9j0RkKnAd1mQiF5SxXrAGIRuENZ77BMckJp406D9LaXtwIY+EzKGZZPKbacKzRTeyI2EQCyf181kZAGyYA9/8HbLTISYRBjwBnW706XPxRBxaRgDGoGVoGR7izS6p07FGv3ynnPXXAMn2rRfwqv3XoyZE/sR1R94kQk4CkCiZPB36Jp9GxWHNHe+bMtgwBz6ZCIX2IJzZe63H4PabGyhxaBkBFoOWoWV4kFeHuRCRJODTcmoKrwOLjTEz7cfbgEvtMfHLlZKSYipznULx/zuf4Jz0ijdUSqmaJKYFPLjR7c1FZLUxJqWi7fzZptCc06f6S6ec6Q1F5C4RSRWR1IyMjEodJDjHG5NjKaWUn2V758dujbii2RjzBvAGWDWFSu0ck2hVt0opaNCMdTf84FYRXeZeTFjub2fGFdMCcTdTv3BBmXFUJtvnPnMeDfLKqEhV5heDB+LQMgIsBi2jjpaReOYyD/BnTWEfp8//mog35rwd8ASE1j99WWh9wq7+G73axLp1C7v6b2eUccLU47vm91Q7DgY84dbuy3ce5vGc4Zww9c6IY3PHB3wWh5YRgDFoGVqGB/kzKSwAxomlN5BdUXtClXS6Ea5/keKoREoQiqMS4foXK9dAU6qM3PCmTC68g9vXtubLTQcq3t+lDGJaAGL9dTOOnRnHuWfGauYV9eXLNo8548gKTWBy4R3cuKwFWw8c83ocWkaAxqBlaBke5LWGZhGZCVwKNMGaNP2vQCiAMeY1u0vqy1jz5J4AbjXGVNiCXNmGZm955bsdPPflNsJDg5hzdx86JTb0ynGO5J5k2H9/ZPfhE1zZMYHXbu5OcJAAUFJimDhrLZ9u2E+zmHDm39eX+Ohwr8ShlKrZ3G1ornGT7ARKUjDG8PCHG/hwdTpxUWHMv68vzRvWr3jHSsgvLObmN1eSuvsoFzaPYfbdvYmoF3LGNjf9bwVr9mTRKTGGWXeduY1SStWE3kc1mojw1LAL6dMmloycAm6fvoqc/EKPlW+M4dG5G0jdfZSmMeG8OT6lzC/78NBg/jcuhRaN67MhPZsHZ6+jpKRmJXqlVODQpFAN9UKCeO3m7rSJa8DWAznc//5aiopLPFL2C4u28/G632hQL5ipE3qQcJbTQrGRYUyb0IOo8BC+3HSQp7/Y6pEYlFJ1jyaFaoqJCGXahB40blCPJb9k8LdPNlPdU3IfrUnnxW+2EyTw8k3d6NA0usJ92sZH8frN3QkJEt74/lfeX7mnWjEopeomTQoe0Cq2AW/c0p16wUG8u2I3U39Mq3JZK389zKNzNwAwZfD5XHZevNv7XtS2CU8NuxCAv3y8ke9/qdyFfkoppUnBQ1KSGvPcyE4APPnZZr7efLDSZezKzOXuGaspLDbc2jeJcX2SKl3GjT1a8LtLz6W4xHDfe2vYdiCn0mUopeouTQoeNKRLc/5wZTuMgYkz17JxX7bb+x7NPcmt034i60QhV3SI5/FrO1Y5joevas+1FzYlp6CI26avIiOnoMplKaXqFk0KHvb7y9syvFtz8gqLuf3tVezPzqtwn4KiYu5+dzVph09wfrNo/jO6q/NahKoIChL+342d6dqyIfuy8rjjnVTyThZXuTylVN2hScHDRIR/Dr+Qnq0bc/BYAbdNT+V4QVG52xtjmDz3Z35KO8I50eG8Nb4HDcKqf52Bo6tqYqP6rN+bxR/maFdVpVTFNCl4QVhIMK/f3J3WTRqwZf8xJs4sv6vqi9/sYN7afUTUC+atCSmcE+O5K5KbuHRV/XzjAZ79cpvHylZK1U6aFLykUYN6TJ3Qg4YRoXy79RBPfrbljG0+XrePFxb9QpDAS2O6cn6zGI/HkZwQxatjra6qry3ZyayftKuqUqp8mhS8qHWTBrxxSwr1goOYviyN6T/ucq5blXaEhz+wup7+5bqODOiQ4LU4Lk5uwpNDrXmOHp+/kR93ZHrtWEqpmk3HPvKBeWvTeXD2+jLXNYoIZe0TV/kkjn9+voXXl/xa5rqOTaMrN8+zUqpG0bGPAsiwrolc2PzMq5IFGHRhU5/F8ejV59Gi8ZmD9oUGC91aNfJZHEqpwKVJwUfeHJdC6V6m9UKCmHRFss9iCAoSZtzei9KdXYNFmDigrc/iUEoFLk0KPpIQU59RPVo4v5BDg4WRKS2Ij/Lt/AetYhswvNupqbBDgoQRfohDKRWYNCn40INXtKNeiPWS+/PX+aMDzzs1UY8xWktQSjlpUvCh+OhwRnZPRAS//jqPjw5ncGerLcMY66aUUqBJwecmDkimR1Jjv/86/9M1HWgUEYoB3lux26+xKKUChyYFH4uPDmfO3X38fg4/Pjqc127uDsD7P+2hoEjHRlJKaVKo03q2bkyHptFkHj/JZxv2+zscpVQA0KRQh4kIEy5qBcD0ZWnVnjFOKVXzaVKo44Z0aU7DiFA2pGezdm+Wv8NRSvmZJoU6Ljw0mNE9WgLw9rI0/wajlPI7TQqKW/q0Ikjgsw37OXQs39/hKKX8SJOConnD+lzV8RyKSgzvrdShtZWqyzQpKAAm9E0C4L2VezhZVPaEQEqp2k+TggKgV+vGnHdOFJnHC1j4s3ZPVaqu0qSgAEf31CQApmmDs1J1liYF5TSkS3Ni6oeyfm8Wa/cc9Xc4Sik/0KSgnOrXC2Z0zxaAdk9Vqq7SpKBOc0tvu3vqz/s5lKPdU5Wqa7yaFERkoIhsE5EdIjK5jPUtReQ7EVkrIhtEZJA341EVS2wUwZUdEygsNryv3VOVqnO8lhREJBh4BbgG6AiMEZGOpTZ7HJhjjOkKjAb+6614lPvG2w3O2j1VqbrHmzWFnsAOY8yvxpiTwCxgSKltDOCY0T4G+M2L8Sg39WkTS/uEKDJyCvh8o3ZPVaou8WZSaA7sdXmcbi9zNQW4WUTSgYXA78sqSETuEpFUEUnNyMjwRqzKhYg4awvTtcFZqTrF3w3NY4DpxphEYBDwroicEZMx5g1jTIoxJiUuLs7nQdZFQ7s2I6Z+KGv3ZLFeR09Vqs7wZlLYB7RweZxoL3N1OzAHwBizHAgHmngxJuWmiHohjOqh3VOVqmu8mRRWAcki0lpE6mE1JC8otc0eYACAiHTASgp6fihAOLqnfrLhNzJyCvwdjlLKB7yWFIwxRcD9wJfAFqxeRptE5O8iMtje7I/AnSKyHpgJTDA6/VfAaNE4ggEdrO6pM3/S7qlK1QVS076DU1JSTGpqqr/DqDOW7cjkpjdXEh8Vxg+PXk69EH83QymlqkJEVhtjUiraTv/D1Vn1OTeWdgmRHMop4ItNB/wdjlLKyzQpqLM6rXvqj7v8G4xSyus0KagKDevanOjwENbsyWJDunZPVao206SgKuTaPVUvZlOqdtOkoNxyS+8kRODT9fvJPK7dU5WqrTQpKLe0jI1gwHnxnCwuYaaOnqpUraVJQbltwkWtAZixcjeFxTp6qlK1kSYF5ba+bWNpGx/JwWMFfLFRu6cqVRtpUlBuc+2equMhKVU7aVJQlTK8a3OiwkNI3X2Ujfuy/R2OUsrDNCmoSmkQFsKNKdo9VanaKsTfAaiaZ8m2QwB8uDqdD1enO5d3bBrNwkn9/BWWUsoDtKagKq13m1ik1LLQYKFbq0Z+iUcp5TmaFFSlTRyQTEjw6WkhWISJA9r6KSKllKdoUlCVFh8dzqiUU5PqBQcJI1JaEB8V7seolFKeoElBVcnEAcmEBFm1hRJjtJagVC2hSUFVSXx0ODd0SwTAGDh0TMdDUqo20KSgquyPV7UjIToM0IvZlKotNCmoKouPDmf2XX0QgY/X/8aR3JP+DkkpVU2aFFS1JDVpwKXt4jhZVMLMn3T0VKVqOk0Kqtom9LVHT12xmyIdPVWpGk2Tgqq2fm2b0KZJA/Zn5/PV5oP+DkcpVQ2aFFS1BQWdGj1Vx0NSqmbTpKA84obuiUSGhfDTriNs/u2Yv8NRSlWRJgXlEZFhIYzobl23oN1Tlaq5NCkojxnXpxUA89ft46h2T1WqRtKkoDymTVwkl7aPo6CohFmr9vo7HKVUFWhSUB7laHB+d3madk9VqgbSpKA8qn9yHK2bNOC37HwWbdHuqUrVNJoUlEcFBYmzbWHaj2n+DUYpVWleTQoiMlBEtonIDhGZXM42N4rIZhHZJCLvezMe5RsjuifSoF4wK3cdYct+7Z6qVE3itaQgIsHAK8A1QEdgjIh0LLVNMvAnoK8x5nzgAW/Fo3wnKjxUu6cqVUN5s6bQE9hhjPnVGHMSmAUMKbXNncArxpijAMaYQ16MR/nQOLvBWbunKlWzuJUUROQjEblWRCqTRJoDrv0S0+1lrtoB7UTkRxFZISIDyzn+XSKSKiKpGRkZlQhB+cu5cZFc0i6O/MISZqdq91Slagp3v+T/C9wEbBeRp0WkvYeOHwIkA5cCY4D/iUjD0hsZY94wxqQYY1Li4uI8dGjlbbc6u6fq6KlK1RRuJQVjzCJjzFigG5AGLBKRZSJyq4iElrPbPqCFy+NEe5mrdGCBMabQGLML+AUrSahaoH+7OJJiI9iXlceiLXpmUKmawO3TQSISC0wA7gDWAv/BShJfl7PLKiBZRFqLSD1gNLCg1DbzsWoJiEgTrNNJv7ofvgpkVvfUJEAbnJWqKdxtU5gHLAUigOuNMYONMbONMb8HIsvaxxhTBNwPfAlsAeYYYzaJyN9FZLC92ZfAYRHZDHwHPGyMOVy9p6QCyYiURCLqBbP818NsPaDdU5UKdGKMqXgjkcuMMd/5IJ4KpaSkmNTUVH+HoSrhiY838s7y3Yzp2ZJ/Dr/Q3+EoVSeJyGpjTEpF27l7+qijawOwiDQSkXurHJ2qUxynkOatTSfrhHZPVSqQuZsU7jTGZDke2NcV3OmdkFRt0zY+kn7JTcgvLGGOdk9VKqC5mxSCRUQcD+yrlet5JyRVG02wu6e+s3w3xSUVn7JUSvmHu0nhC2C2iAwQkQHATHuZUm65rH08rWIjSD+axzc6eqpSAcvdpPAoVu+g39m3b4BHvBWUqn2CgoRbelujp07X7qlKBSx3L14rMca8aowZYd9eN8YUezs4VbuMTGlBRL1glu08zC8Hc/wdjlKqDO5ep5AsIh/aQ1z/6rh5OzhVu8TUD2V4N2v4K60tKBWY3D19NA14FSgCLgPeAWZ4KyhVe413dE9ds4/sE4X+DUYpdQZ3k0J9Y8w3WBe77TbGTAGu9V5YqrZKToji4rZNyCss1u6pSgUgd5NCgT1s9nYRuV9EhlHO8BZKVcTZPXVFmnZPVSrAuJsUJmGNezQR6A7cDIz3VlCqdrvsvHhaNK7P3iN5fLtVR09VKpCEVLSBfaHaKGPMQ8Bx4FavR6VqteAg4WSRNb/Cne+cPo5Vx6bRLJzUzx9hKaVwo6Zgdz292AexqDrkknZnTpYUGix0a9XID9EopRwqrCnY1orIAuADINex0BjzkVeiUrXew1e156PV6RS7NCkEizBxQFv/BaWUcjsphAOHgctdlhlAk4KqkvjocK7smMAXm6whL0KDhREpLYiPCvdzZErVbW4lBWOMtiMoj/vToA7OpBCktQSlAoJbSUFEpmHVDE5jjLnN4xGpOqNVbAMiw4I5XlDMlR0TtJagVABw9/TRpy73w4FhwG+eD0fVNd1bNmLJ9kx6tW7s71CUUrh/+miu62MRmQn84JWIVJ3SqUVDlmzP5MCxfH+HopTC/YvXSksG4j0ZiKqbkhOiAPjl4HE/R6KUAvfbFHI4vU3hANYcC0pVS3tnUtChtJUKBO6ePorydiCqbmrdpAEhQcKeIyc4cbKIiHruNnMppbzB3fkUholIjMvjhiIy1HthqbqiXkgQbeIaYAzsOKSnkJTyN3fbFP5qjMl2PDDGZAF/9U5Iqq5pZ59C2nZATyEp5W/uJoWyttN6vvIIR1LYrjUFpfzO3aSQKiLPi8i59u15YLU3A1N1h9YUlAoc7iaF3wMngdnALCAfuM9bQam6pf052gNJqUDhbu+jXGCyl2NRdVTLxhGEhQSxPzuf7LxCYuqH+jskpeosd3sffS0iDV0eNxKRL70XlqpLgoOE5ARrdtftWltQyq/cPX3UxO5xBIAx5ih6RbPyoHbxemWzUoHA3aRQIiItHQ9EJIkyRk1VqqraabuCUgHB3aTwZ+AHEXlXRGYAS4A/VbSTiAwUkW0iskNEym2TEJEbRMSISIqb8ahapr32QFIqILiVFIwxXwApwDZgJvBHIO9s+4hIMPAKcA3QERgjIh3L2C4KmASsrFTkqlbRmoJSgcHdhuY7gG+wksFDwLvAlAp26wnsMMb8aow5idWVdUgZ2/0f8AxWN1dVRzWLCScyLITDuSfJPF7g73CUqrPcPX00CegB7DbGXAZ0BbLOvgvNgb0uj9PtZU4i0g1oYYz57GwFichdIpIqIqkZGRluhqxqEpFTPZC0tqCU/7ibFPKNMfkAIhJmjNkKtK/OgUUkCHgeq/ZxVsaYN4wxKcaYlLi4uOocVgUw5zDa2q6glN+4O35Run2dwnzgaxE5CuyuYJ99QAuXx4n2Moco4AJgsYgAnAMsEJHBxphUN+NStYhzuAvtlqqU37h7RfMw++4UEfkOiAG+qGC3VUCyiLTGSgajgZtcyswGmjgei8hi4CFNCHWXDnehlP9VeqRTY8wSN7crEpH7gS+BYGCqMWaTiPwdSDXGLKjssVXt1s5lFjZjDHYNUinlQ14d/toYsxBYWGrZE+Vse6k3Y1GBr0lkPRpFhHL0RCEHjuXTNKa+v0NSqs5xt6FZKa8TER1GWyk/06SgAoq2KyjlX5oUVEA5VVPQHkhK+YMmBRVQHDWF7Ye0pqCUP2hSUAHl1BDaOZSU6EC8SvmaJgUVUGIiQkmIDiO/sIS9R0/4Oxyl6hxNCirgaA8kpfxHk4IKOO0TtAeSUv6iSUEFnFNzK2gPJKV8TZOCCjjttKaglN9oUlABJznemldhZ8ZxCotL/ByNUnWLJgUVcBqEhdCicX0Kiw1pmbn+DkepOkWTggpI7Z1zK+gpJKV8SZOCCkin2hW0sVkpX9KkoAJSO52aUym/0KSgApL2QFLKPzQpqIDUJq4BwUFC2uFc8guL/R2OUnWGJgUVkMJDg0mKjaDEwI5D2q6glK9oUlABS4fRVsr3NCmogJUcrxPuKOVrmhRUwNKpOZXyPU0KKmDpENpK+Z4mBRWwkmIjqBccxL6sPI4XFPk7HKXqBE0KKmCFBAdxrj043nY9haSUT2hSUAGtXYKVFLRdQSnf0KSgAtqpdgXtgaSUL2hSUAFNp+ZUyrc0KaiApt1SlfItTQoqoDVvWJ+IesEcyingaO5Jf4ejVK2nSUEFtKAgcU7PqbUFpbwvxJuFi8hA4D9AMPCmMebpUuv/ANwBFAEZwG3GmN3ejEnVPO0Solifns0vB3Po1SbWo2UP+s9SNu8/dsbyjk2jWTipn9f3VyrQeK2mICLBwCvANUBHYIyIdCy12VogxRjTCfgQeNZb8aiay9Gu4I2pObu1bEhosJy2LDRY6NaqkU/2VyrQeLOm0BPYYYz5FUBEZgFDgM2ODYwx37lsvwK42YvxqBrKm1NzThyQzAer0wHjXFZSYsg+cZKHP1hf4f55hcUUl5jTlgWLMHFAW0+HqpRPeDMpNAf2ujxOB3qdZfvbgc/LWiEidwF3AbRs2dJT8akawrUHkjEGEalgD/fFR4fTKjbitIRTbOCTDfurVF6QwIiUFsRHhXsqRKV8yqttCu4SkZuBFKB/WeuNMW8AbwCkpKSYsrZRtVd8VBjR4SFknSgkI6eA+GjPfeEeyy9k75ETzsehwcLkgecRFR7qfhl5hTz9xVaKSgwlBsb1aeWx+JTyNW8mhX1AC5fHifay04jIFcCfgf7GmAIvxqNqKBGh/TlRrEo7yraDOR5NCh+mppNXWEJ8VBgZxwsY1aMlt/drU+ly0g7nMmPlHgAWbTnoPOWlVE3jzS6pq4BkEWktIvWA0cAC1w1EpCvwOjDYGHPIi7GoGs4bw2iXlBjeWZ4GwB+vbEePpMZVbguYOCDZOU7TjOW7KSou8VCUSvmW15KCMaYIuB/4EtgCzDHGbBKRv4vIYHuz54BI4AMRWSciC8opTtVxzqk5PdjYvOSXDNIOn6B5w/rc0D2ROXf3qXJbQHx0OF9MuoQ2TRrwW3Y+X28+6LE4lfIlr7YpGGMWAgtLLXvC5f4V3jy+qj2cNQUPdkudviwNgJt7tyIkuPq/j4KChHF9WjHlk81MX5bGNRc2rXaZSvmaXtGsagRHUth+MIeSkur3NdiZcZwlv2QQFhLE6B4tKt7BTTd0TyQyLISVu46wpYyL2pQKdJoUVI3QuEE9mkSGkXuymH1ZedUu793l1oXzQ7s0p1GDetUuzyEqPJQR3RMBeNuuiShVk2hSUDVG+3M8MwZSTn4hH6Ral9CMvyipumGdwdEldd7afTqIn6pxNCmoGsNTVzbPXZ1O7slierZuTMdm0Z4I7TRt4iLp3y6OgqISZqfurXgHpQKIJgVVY3hiwp2SEsPb9qmjCV6oJTg4yn5Xu6eqGkaTgqoxkj1wrcL32zPYlZlL05hwruqY4KnQztC/XRxJsRHsy8pj0Ra9BEfVHAExzEV1FRYWkp6eTn5+vr9DUW4KDw8nMTGR0FD3h5NwXBy2I+M4xSWG4KDKj4Hk6W6o5bG6pybx9083M33ZLgZecI7XjqWUJ9WKpJCenk5UVBRJSUkeHSxNeYcxhsOHD5Oenk7r1q3d3i8qPJTmDeuzLyuP3YdzaRMXWanj7srMZfG2DOqFBDGmp/cHVhyRksj/+2obK349wtYDxzjvHM+3XyjlabXi9FF+fj6xsbGaEGoIESE2NrZKNTtHbaEq7QqOIS23j3VgAAAbgUlEQVSGdG5GYw92Qy1PdHgoNzi7p+rcUapmqBVJAdCEUMNU9f1q55hw50DleiAdLyjig9R0wDvdUMszro91rHlr08k6od1TVeCrNUlB1Q3t4qvWA+mjNekcLyiiR1IjLmge443QytQ2PpJ+yU3ILyxhjnZPVTVArWhTqAxvzal70UUXsWzZskrvl5SURGpqKk2aNCl3m6eeeorHHnus2scqbevWrYwePRoR4cMPP+Tcc8+tdpmlzZ8/n3bt2tGxY+mZWKumKlNzlpQYZwPzhIvcb8PwlFv7JrF0eybvLN/N7Re3qVIDuS/ofNMK6mBNwVtz6nriS7o8Tz31lFeONX/+fEaMGMHatWvdSgjGGEpKKtfnfv78+WzevLniDd3UNj4SEUjLzKWgqNitfX7YkcmvGbmcEx3OVed7rxtqeS5tF0+r2AjSj+bxzZbAHT1V55tWUAtrCkmTP6v0PoXFhhkrdjNjRfmNgWlPX3vWMiIjIzl+/Dj79+9n1KhRHDt2jKKiIl599VX69evHzJkzeeqppzDGcO211/LMM8+cUcbQoUPZu3cv+fn5TJo0ibvuuovJkyeTl5dHly5dOP/883nvvfecxzLG8Mgjj/D5558jIjz++OOMGjWKxYsXM2XKFJo0acLGjRvp3r07M2bMOO08/sKFC/n3v/9NcHAw33zzDd999x3PP/88U6dOBeCOO+7ggQceIC0tjauvvppevXqxevVqFi5cyLZt2/jrX/9KQUEB5557LtOmTSMyMpLJkyezYMECQkJCuOqqqxg+fDgLFixgyZIlPPnkk8ydO7fatZHw0GCSYhuwKzOXXZm5bvXocdQSbunTilAvdkMtj6N76v99ao2eetX5gdk9taz5qnW+6bqn1iUFf3v//fe5+uqr+fOf/0xxcTEnTpzgt99+49FHH2X16tU0atSIq666ivnz5zN06NDT9p06dSqNGzcmLy+PHj16cMMNN/D000/z8ssvs27dujOO9dFHH7Fu3TrWr19PZmYmPXr04JJLLgFg7dq1bNq0iWbNmtG3b19+/PFHLr74Yue+gwYN4p577iEyMpKHHnqI1atXM23aNFauXIkxhl69etG/f38aNWrE9u3befvtt+nduzeZmZk8+eSTLFq0iAYNGvDMM8/w/PPPc9999zFv3jy2bt2KiJCVlUXDhg0ZPHgw1113HSNGjPDYa9wuIZJdmblsO5BTYVJIy8zlu22HqOfh0VAra6TdPXXZzsP8cjAnIGdmi48Op1NiDKvSjjqXDevWXOebrmNqXVKo6Bc9wKFj+fR79jsKikoIDwni+0cv89gHv0ePHtx2220UFhYydOhQunTpwrfffsull15KXFwcAGPHjuX7778/Iym8+OKLzJs3D4C9e/eyfft2YmNjyz3WDz/8wJgxYwgODiYhIYH+/fuzatUqoqOj6dmzJ4mJVnfILl26kJaWdlpSKKusYcOG0aBBAwCGDx/O0qVLGTx4MK1ataJ3794ArFixgs2bN9O3b18ATp48SZ8+fYiJiSE8PJzbb7+d6667juuuu66Kr2DF2iVE8eWmg241Nr+zfDfGwODOzYiNDPNaTBWJDg/lhm6JvLtiN9OXpfHUsAv9Fkt5vtt6iFSXhADW/0pJiSEoQNtBlOfVuTYFsH4RjeyeiAiMSGnh0V9Cl1xyCd9//z3NmzdnwoQJvPPOO27tt3jxYhYtWsTy5ctZv349Xbt2rdYV2mFhp74Ag4ODKSoqqnJZjkQBVrvClVdeybp161i3bh2bN2/mrbfeIiQkhJ9++okRI0bw6aefMnDgwCofryKnpuY8e7fU3IIi52io3hznyF3jL7JHT12zj+wThX6O5nSbfzvG/e+vwQAXNotGgNAg4ZutGTz31TZ/h6d8qE4mBbDOn1ZnTt7y7N69m4SEBO68807uuOMO1qxZQ8+ePVmyZAmZmZkUFxczc+ZM+vfvf9p+2dnZNGrUiIiICLZu3cqKFSuc60JDQyksPPNLpF+/fsyePZvi4mIyMjL4/vvv6dmzZ5Xi7tevH/Pnz+fEiRPk5uYyb948+vU7s8dJ7969+fHHH9mxYwcAubm5/PLLLxw/fpzs7GwGDRrECy+8wPr16wGIiooiJ8dzs6WBy9Sch85e7kdr0skpKCKllW+7oZanbXwU/ZKbkFdYHFDdUw8ey+f2t1eRe7KYwZ2b8eb4FHq0bszzozoTHCS8ungnc1YFTrzKu+psUoiPDq/WnLzlWbx4MZ07d6Zr167Mnj2bSZMm0bRpU55++mkuu+wyOnfuTPfu3RkyZMhp+w0cOJCioiI6dOjA5MmTnadrAO666y46derE2LFjT9tn2LBhdOrUic6dO3P55Zfz7LPPcs45VWvE7NatGxMmTKBnz5706tWLO+64g65du56xXVxcHNOnT2fMmDF06tSJPn36sHXrVnJycrjuuuvo1KkTF198Mc8//zwAo0eP5rnnnqNr167s3LmzSrGVlhTbgNBgYc+RE5w4WXYNyJhT3VB9ebFaRcbbF7O9syKNYg/MIFddJ04Wcfvbq9ifnU9Kq0Y8O6ITCTH1mXN3H67v3Jwnh14AwGPzfubHHZl+jlb5ghjj/w9mZaSkpJjU1NTTlm3ZsoUOHTr4KSJVVdV5365+4Xu2Hcxhwf196ZTY8Iz1S7dncMtbP5EQHcYPj17ul15HZSkuMVz2r8XsOXKC/41L4UovjtTqTiz3zFjN15sP0rJxBPPuvajMdpd/LtzC69//SlR4CPPuvYi28YHXSK4qJiKrjTEpFW0XGP8pSlVSsj0GUnnDaDumwry5l3+6oZYnOEicM7P5e7rOfy7cwtebDxIdHsLUCT3KbYh/dOB5DDz/HHLyi7h1+ioyjxf4OFLlS4Hz36JUJZxtwp09h0/wzdZD1AsOYkwv74+GWlkjU1pQPzSYH3Zksr2aU4tW1YwVu3nzh12EBAmv3dKdtvHljzgbFCS8MKoLnRNj2Hskj7veSSW/0L0LB1XNo0lB1UiOgfHKmprzneVpGAPXdW5KEz92Qy1PTP1QhndrDsDby9N8fvwlv2Tw1wWbAPjn8Au56Nzyh1hxqF8vmP+NT6F5w/qs2ZPFQx+spyQA2kSU52lSUDVSeTWF3IIi57zIt/phnCN3ObrIzl29j+w833VP3XrgGPe9t4biEsN9l53LyBT3L+iLjwpn6oQeRIaF8OmG/Tz/9S9ejFT5iyYFVSO1aBxBeGgQ+7PzT/tSnbd2Hzn5RXRr2ZALE/3fDbU8yQlR9G0bS15hsfNaCm87lJPP7dNTOV5QxHWdmvLHK9tXuoz250TxythuBAcJL3+3w2exK9/RpKBqpOAgcZ4Hd5yXN8Y4G28n9A3cWoKDY8TWd5bv9nr31LyTxdz5dir7svLo2rIh/xrZucpXKfdvF8ffBp8PWF1Vl+887MlQlZ/VzaSwYQ68cAFMaWj93TDHo8VPmTKFf/3rXx4tc9CgQWRlZZGVlcV///vfSu//8MMPc/755/Pwww97NC6HqsZVHc4rm+2ksGznYbYfOk58VBjX1IA5kS8/L57ERvXZc+QEi7cd8tpxSkoMD85ex/r0bFo0rs//xqUQHhpcrTJv7t2KOy5uTWGx1a11Z0blJj1SgavuJYUNc+CTiZC9FzDW308mejwxeNrChQtp2LBhlb9833jjDTZs2MBzzz3n1vaVHRbDH0nB0a6w3W5snvZjGmB9YQVSN9TyBAeJ82K26V7snvrMF1v5YtMBosJDmDahh8ca3/80qANXdkwgO6+Q26av4kiuzixXGwT+f05lTYk5++2jO6Ew7/R9CvOs5WfbrwL/+Mc/aNeuHRdffDHbtp0aK2bnzp0MHDiQ7t27069fP7Zu3QrAhAkTmDhxIhdddBFt2rThww8/BGD//v1ccskldOnShQsuuIClS5cC1mQ8mZmZTJ48mZ07d9KlSxcefvhhxo0bx/z5853HGzt2LB9//PFpsQ0ePJjjx4/TvXt3Zs+eTVpaGpdffjmdOnViwIAB7NmzxxnTPffcQ69evXjkkUfIzc3ltttuo2fPnnTt2tVZ7qZNm+jZsyddunShU6dObN++/Yy4fOHU1Jw57D1ygm+2HrS6ofYMvG6o5bnR7p66dHsmOyoYtqMqZv60h9e//9Xqenpzd49eeBYcJPxndBcubB7D7sMntKtqLVHrRkn1h9WrVzNr1izWrVtHUVER3bp1o3v37oA1RMVrr71GcnIyK1eu5N577+Xbb78FrATwww8/sHXrVgYPHsyIESPKHHrb1dNPP83GjRudQ2kvWbKEF154gaFDh5Kdnc2yZct4++23T9tnwYIFREZGOve5/vrrGT9+POPHj2fq1KlMnDjRmVjS09NZtmwZwcHBPPbYY1x++eVMnTqVrKwsevbsyRVXXMFrr73GpEmTGDt2LCdPnqS4uPiMuHzBtQeSsxtqp6bERQVeN9TyxESEMqxbc95fuYe3l+3m/+xhJTxh6fYMHp+/EYB/DLuAvm0r7npaWRH1QnhzfApDX/mR1N1HeXTuBv49qovOmV6D1b6kMCX77OtfuMA+dVRKTAt4cGOVDrl06VKGDRtGREQEYP0yBzh+/DjLli1j5MiRzm0LCk5dDTp06FCCgoLo2LEjBw9aM3KVNfT22fTv3597772XjIwM5s6dyw033EBIyNnf1uXLl/PRRx8BcMstt/DII484140cOZLgYOt881dffcWCBQuc7SP5+fns2bOHPn368I9//IP09HSGDx9OcnKyW6+TpzWNCScqLITDuSd5b6VV2wmkcY7cNb5PEu+v3MPcNek8PLA90eGh1S7zl4M53DvD6np6T/9zGdXDe7WnhGirq+qIV5fx8brfaBXbgD9c2c5rx1Pe5dWkICIDgf8AwcCbxpinS60PA94BugOHgVHGmDRvxsSAJ6w2BNdTSKH1reUeVlJSQsOGDcv99ew6vLVjDCrH0NufffYZEyZM4A9/+APjxo0763HGjRvHjBkzmDVrFtOmTatWzKWHyZ47dy7t25/edbFDhw706tWLzz77jEGDBvH666/Tpk2bah23Kq598QdyCqy2jxMnrdMWQ175scbNKfzgbOvzceJkMZ2mfOVcXpnnUd78ytHhITxydeW7nlZWh6bRvHxTN26dvooXv9nOi99sP229J56LllH1MirDa20KIhIMvAJcA3QExohI6dnbbweOGmPaAi8AZ85R6WmdboTrX7RqBoj19/oXreVVdMkllzB//nzy8vLIycnhk08+ASA6OprWrVvzwQcfANaXrGNI6fKUNfS2q7KGop4wYQL//ve/AejYsfRLfKaLLrqIWbNmAfDee++VOUQ2wNVXX81LL73kTFhr164F4Ndff6VNmzZMnDiRIUOGsGHDBq8MkV2Rbi0bUrpXZU2cU7hby4YEB505N/L5zaI4fLzArVvHZlFnzK8McG2npj6bIOey8+LpkXTma++J56Jl+G7ubG/WFHoCO4wxvwKIyCxgCOA6i/sQYIp9/0PgZRER4+2hWzvdWK0kUFq3bt0YNWoUnTt3Jj4+nh49ejjXvffee/zud7/jySefpLCwkNGjR9O5c+dyy1q8eDHPPfccoaGhREZGnjFJT2xsLH379uWCCy7gmmuu4bnnniMhIYEOHTqcMZNbeV566SVuvfVWnnvuOeLi4sqtXfzlL3/hgQceoFOnTpSUlNC6dWs+/fRT5syZw7vvvktoaCjnnHMOjz32GI0bNz4jLm+bOCCZWav2UmJq9pzCjrmRXa9VKCw2fLB6Hx+s3lflcsNCgnjQx6dxXrmpG33++S3FxrPPRcs4nTc/514bOltERgADjTF32I9vAXoZY+532WajvU26/XinvU1mqbLuAu4CaNmyZffdu3efdqy6PnT2iRMnuPDCC1mzZg0xMYF7FW9pnnjfHpy9lnlrfwOsX0+jerR0zgFQkzw+72dmrtpDcYn1OCwkiAZhlfvNlltQREGRVUBIkDC6p39eiz/P+5mZP+3BkeOq+1y0jNPLqOrn3N2hs2tEQ7Mx5g3gDbDmU/BzOAFl0aJF3H777Tz44IM1KiF4yp+u6cDCnw9QUFRSI2sJDqdqC1WfN9x17vGQIP+9FpMGJPPh6vRqzYHuiXnUa2sZ3v6ce/M6hX2A62hbifayMrcRkRAgBqvBWbnpiiuuYPfu3TzwwAP+DsUvvDnfti954nkEymsRKM9Fy6gab9YUVgHJItIa68t/NHBTqW0WAOOB5cAI4NuqticYY7RvdA3iydOWEwck88uh4zW2luDgiecRKK9FoDwXLaPyvDodp4gMAv6N1SV1qjHmHyLydyDVGLNARMKBd4GuwBFgtKNhujxlTce5a9cuoqKiiI2N1cRQAxhjOHz4MDk5ObRuHfgD1ylVG7jbplAr5mguLCwkPT2d/Px8P0WlKis8PJzExERCQ6t/oZZSqmK1qqG5IqGhofqLUymlPKD2DYinlFKqyjQpKKWUctKkoJRSyqnGNTSLSAawu8INy9YEyKxwK+/TOE6ncQRWDKBxlFYb4mhljImraKMalxSqQ0RS3Wl91zg0jrocg8ZRt+PQ00dKKaWcNCkopZRyqmtJ4Q1/B2DTOE6ncZwSCDGAxlFanYmjTrUpKKWUOru6VlNQSil1FpoUlFJKOdWqpCAiU0XkkD2jm2NZYxH5WkS2238b2ctFRF4UkR0iskFEunkohhYi8p2IbBaRTSIyyU9xhIvITyKy3o7jb/by1iKy0j7ebBGpZy8Psx/vsNcneSIOl3iCRWStiHzqrzhEJE1EfhaRdSKSai/z6ftil91QRD4Uka0iskVE+vjh89Hefh0ct2Mi8oAf4njQ/nxuFJGZ9ufWH5+NSXYMm0TkAXuZ118L8dB3loiMt7ffLiLjq/5KYA1jXFtuwCVAN2Cjy7Jngcn2/cnAM/b9QcDngAC9gZUeiqEp0M2+HwX8AnT0QxwCRNr3Q4GVdvlzsIYoB3gN+J19/17gNfv+aGC2h9+bPwDvA5/aj30eB5AGNCm1zKfvi13228Ad9v16QEN/xOESTzBwAGjlyziA5sAuoL7LZ2KCrz8bwAXARiACa5DQRUBbX7wWeOA7C2gM/Gr/bWTfb1TlmDz9AfP3DUgq9QJvA5ra95sC2+z7rwNjytrOw/F8DFzpzzjsD/saoBfW1ZAh9vI+wJf2/S+BPvb9EHs78dDxE4FvgMuBT+0PtT/iSOPMpODT9wVrdsFdpZ+Tnz8fVwE/+joOrKSw1/4yC7E/G1f7+rMBjATecnn8F+ARX70WVPM7CxgDvO6y/LTtKnurVaePypFgjNlv3z8AJNj3HR9Ih3R7mcfY1duuWL/SfR6HfcpmHXAI+BrYCWQZY4rKOJYzDnt9NhDriTiwJlp6BHDMXh7rpzgM8JWIrBaRu+xlvn5fWgMZwDT7dNqbItLAD3G4Gg3MtO/7LA5jzD7gX8AeYD/We70a3382NgL9RCRWRCKwfpG3wH/vSWWP69F46kJScDJWGvVJH1wRiQTmAg8YY475Iw5jTLExpgvWL/WewHnePmZpInIdcMgYs9rXxy7DxcaYbsA1wH0iconrSh+9LyFYpwteNcZ0BXKxThH4Og4A7PP1g4EPSq/zdhz2ufIhWImyGdAAGOit45XHGLMFeAb4CvgCWAcUl9rGZ++Jv49bF5LCQRFpCmD/PWQv34f1a8Ah0V5WbSISipUQ3jPGfOSvOByMMVnAd1hV8YYi4phcyfVYzjjs9THAYQ8cvi8wWETSgFlYp5D+44c4HL9MMcYcAuZhJUpfvy/pQLoxZqX9+EOsJOGvz8c1wBpjzEH7sS/juALYZYzJMMYUAh9hfV788dl4yxjT3RhzCXAUqy3QX+9JZY/r0XjqQlJYADha48djneN3LB9nt+j3BrJdqmxVJiICvAVsMcY878c44kSkoX2/Pla7xhas5DCinDgc8Y0AvrV/pVSLMeZPxphEY0wS1mmKb40xY30dh4g0EJEox32s8+gb8fH7Yow5AOwVkfb2ogHAZl/H4WIMp04dOY7nqzj2AL1FJML+v3G8Fj79bACISLz9tyUwHKtThL/ek8oe90vgKhFpZNe+rrKXVU11G2kC6Yb14d4PFGL9Irsd65zjN8B2rF4Fje1tBXgF6zz7z0CKh2K4GKu6twGrGroO6xylr+PoBKy149gIPGEvbwP8BOzAOmUQZi8Ptx/vsNe38cL7cymneh/5NA77eOvt2ybgz/Zyn74vdtldgFT7vZmP1WPEH3E0wPqlHeOyzNef078BW+3P6LtAmD8+o8BSrIS0Hhjgq9cCD31nAbfZr8sO4NbqvBY6zIVSSimnunD6SCmllJs0KSillHLSpKCUUspJk4JSSiknTQpKKaWcNCkoVQn29R8r7SEq+vkxjiki8pC/jq9qr5CKN1FKuRgA/GyMucPfgSjlDVpTUDWKiCSJNQfB/8Qa+/4r+4ptRGSxiKTY95vYQ2sgIhNEZL49Nn2aiNwvIn+wf+2vEJHG5RznW3vc+m9EpKWIdMEa1niIWHMQ1C+1z9NizaOxQUT+ZS+73qVmsUhEEuzlU0TkbRFZKiK7RWS4iDwr1nwPX9hDpTjmgHAs/0lE2pYR67n2Pqvt8s6zl48Ua46A9SLyvQffBlWLaVJQNVEy8Iox5nwgC7jBjX0uwBq+oAfwD+CEsQakWw6MK2P7l4C3jTGdgPeAF40x64AnsMbx72KMyXNsLCKxwDDgfHufJ+1VPwC97WPNwhot1uFcrLGgBgMzgO+MMRcCecC1Lttl28tfxhpxtrQ3gN8bY7oDDwH/tZc/AVxtjOlsH0OpCunpI1UT7bK/oMEaajnJjX2+M8bkADkikg18Yi//GWtIkNL6YCURsIZfeLaC8rOBfOAtsWaX+9RengjMtgc2q4c1l4LD58aYQhH5GWuimy9cYnJ9TjNd/r7gelCxRuO9CPjAGj4IsIaKAPgRmC4ic7AGm1OqQlpTUDVRgcv9Yk79uCni1Gc6/Cz7lLg8LsEDP46MNb5/T6xRT6/j1Bf8S8DL9i/9u0vFVWDvWwIUmlNjzpSOyZRzH6znm2XXXBy3Dna59wCPY42gudquzSh1VpoUVG2SBnS37484y3buWIY1qivAWKwB08pl/2KPMcYsBB4EOturYjg1jPH4svZ1wyiXv8tdVxhrro5dIjLSjkNEpLN9/1xjzEpjzBNYE/u4Dq+sVJn09JGqTf4FzBFrVrXPqlnW77FmR3sY6wv11gq2jwI+FpFwrNEs/2Avn4J1auco8C3WhDKV1UhENmDVLMaUsX4s8KqIPI41H/csrNE+nxORZDueb+xlSp2VjpKqVACze1ClGGMy/R2Lqhv09JFSSiknrSkopZRy0pqCUkopJ00KSimlnDQpKKWUctKkoJRSykmTglJKKaf/D6FAT9KOBOkJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xlabel = np.arange(100, 1050, step=50)\n",
    "plt.title('Affects of masking')\n",
    "plt.plot(xlabel, res_iso / 30, label = \"isolation forest\", linewidth = 2, marker = \"v\")\n",
    "plt.plot(xlabel, res / 30, label = \"density forest\", linewidth = 2, marker = \"o\")\n",
    "plt.xticks(np.arange(100,1050, step = 100))\n",
    "plt.legend()\n",
    "plt.ylabel('accuracy')\n",
    "plt.xlabel('num of samples')\n",
    "plt.savefig(\"masking.pdf\", bbox_inches = 'tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This image is the average accuricy when there are 970 random vectors of dimension 10, and 30 points of all zeros. The error rate is set to 5%.  We check what percentage of the 30 zeros were caught by the algorithm. We see the isolation forest is very sensitive to the number of samples per tree (we have 20 trees). Nothing special about these parameters, the experiment results are pretty robust, in particular, the number of trees could be increased significantly without changing the outcome."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
